Fast and stable method for simulating quantum electron dynamics 
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A fast and stable method is formulated to compute the time evolution of a wavefunction by 
numerically solving the time-dependent Schrodinger equation. This method is a real space/real time 
evolution method implemented by several computational techniques such as Suzuki's exponential 
product, Cayley's form, the finite differential method and an operator named adhesive operator. 
This method conserves the norm of the wavefunction, manages periodic conditions and adaptive 
mesh refinement technique, and is suitable for vector- and parallel-type supercomputers. Applying 
this method to some simple electron dynamics, we confirmed the efficiency and accuracy of the 
method for simulating fast time-dependent quantum phenomena. 
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I. INTRODUCTION 



There are many computational method of solving the 
TD-Schrodinger equation numerically. Conventionally, 
a wavefunction has been represented as a linear com- 
bination of plane waves or atomic orbitals. However, 
these representations entail high computational cost to 
calculate the matrix elements for these bases. The plane 
wave bases set is not suitable for localized orbitals, and 
the atomic orbital bases set is not suitable for spreading 
waves. Moreover, they are not suitable for paralleliza- 
tion, since the calculation of matrix elements requires 
massive data transmission among processors. 

To overcome those problems, some numerical methods 
adopted real space representation jl]-f|]. In those meth- 
ods, a wavefunction is dcscritized by grid points in real 
space, and with them some dynamic electron phenomena 
were simulated successfully 

Among these real space methods, a method called 
Cayley's form or Crank-Nicholson scheme is known to 
be especially useful for one-dimensional closed systems 
because this method conserves the norm of the wave- 
function exactly and the simulation is rather stable and 
accurate even in a long time slice. These characteris- 
tics are very attractive for simulations over a long time 
span. Unfortunately, this method is not suitable for two- 
or three-dimensional systems. This problem is fatal for 
physically meaningful systems. Though there are many 
other computational methods that can manage two- or 
three-dimensional systems, these methods also have dis- 
advantages. 

In the present work, we have overcome the problems 
associated with Cayley's form and have formulated a 
new computational method which is more efficient, more 
adaptable and more attractive than any other ordinary 
methods. 

In our method, all computations are performed in real 
space so there is no need of using Fourier transform. The 
time evolution operator in our method is exactly unitary 
by using Cayley's form and Suzuki's exponential product 



so that the norm of the wavefunction is conserved during 
the time evolution. Stability and accuracy are improved 
by Cayley's form so we can use a longer time slice than 
those of the other methods. Cayley's form is a kind of 
implicit methods, this is the key to the stability, but im- 
plicit methods are not suitable for periodic conditions 
and parallelization. We have avoided these problems by 
introducing an operator named adhesive operator. This 
adhesive operator is also useful for adaptive mesh refine- 
ment technique. 

Our method inherits many advantages from many ordi- 
nary methods, and yet more improved in many aspects. 
With these advantages, this method will be useful for 
simulating large-scale and long-term quantum electron 
dynamics from first principles. 

In section II, we formulate the new method step by 
step. In section III, we apply it to some simulations of 
electron dynamics and demonstrate its efficiency. In sec- 
tion IV, we draw some conclusions. 



II. FORMULATION 



In this section, we formulate the new method step 
by step from the simplest case to complicated cases. 
Throughout this paper, we use the atomic units h = 
1, m = 1, e = 1. 



A. One-dimensional closed free system 



For the first step, we consider a one-dimensional closed 
system where an electron moves freely but never leaks 
out of the system. The TD-Schrodinger equation of this 
system is simply given as 



. dip(x, t) 

1 dt 



ip(x,t) 



(1) 
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The solution of Eq. 
ponential operator as 



Gj) is analytically given by an ex- decomposed into the LU form as 



tp(x, t + At) = exp 



iAt- 



ip(x,t) 



(2) 



where At is a small time slice. By using Eq. (Q) repeat- 
edly, the time evolution of the wavefunction is obtained. 

An approximation is utilized to make a concrete form 
of the exponential operator. We have to be careful not to 
destroy the unitarity of the time evolution operator, oth- 
erwise the wavefunction rapidly diverges. We adopted 
Cayley's form because it is unconditionally stable and 
accurate enough. Cayley's form is a fractional approxi- 
mation of the exponential operator given by 



exp 



[At 



<i 2 



l + iAtd 2 /4 
l-iAtflg/4 



(3) 



It is second-order accurate in time. By substituting 
Eq. (JsJ) for Eq. (||) and moving the denominator onto the 
left-hand side, the following basic equation is obtained: 



.At<% 
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ip(x,t + At) = 



.Atflg 
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Here bi and Ui 
bi(t) = 



(8) 



are auxiliary vectors defined as below 

V»i-i(i)+BV<(*)+V'i+i(*). (9) 
l/(A-Ui-i), u = (10) 



The auxiliary vector Ui is determined in advance, and it 
is treated as a constant vector in Eq. ( J10[) . 26N float- 
ing operations are heeded to solve Eq. (|10[); here N is 
the number of the grid points in the system, about twice 
that of the Euler method. Unlike the Euler method, it ex- 
actly conserves the norm because the matrices in Eq. (|^) 
are unitary. Moreover, the expected energy is conserved 
because the time evolution operator commutes with the 
Hamiltonian in this case. 



This is identical with the well-known Crank-Nicholson 
scheme. The wavefunction is descritized by grid points 
in real space as 

ipi{t) = ip(xi,t) ; Xl = iAx, i = 0,--- ,N-1 (5) 

where Ax is the span of the grid points. We approximate 
the spatial differential operator by the finite difference 
method (FDM). Then Eq. (||) becomes a simultaneous 
linear equation for the vector quantity ipi (t + At) . For 
example, in a system with six grid points, Eq. (Q) is ap- 
proximated in the following way: 



A -1 

-1 A -1 
0-14-1 

0-14 



ipi(t + At) 
Mt + At) 
Mt + At) 
^(t + At) 

B 1 
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1 B 









Mt) 











(6) 



In the above, 



Ax 2 

A = -4& + 2 . B 
At 



..Ax 2 „ 

-4i— 2 

At 



(7) 



and a- n d ^5 are fixed at zero due to the boundary 
condition. 

It is easy to solve this simultaneous linear equation be- 
cause the matrix appearing on the left-hand side is easily 



B. Three-dimensional closed free system 

It is easy to extend this technique to a three- 
dimensional system. The formal solution of the TD- 
Schrodinger equation in a three-dimensional system is 
given by an exponential of the sum of three second dif- 
ferential operators as 



tp(r, t + At) = exp 



■a ( d l d l d l 
At + + -2- 

V 2 2 2 



V(r,t) . (11) 



These differential operators in Eq. ( |TT| ) are commutable 
among each other, so the exponential operator is exactly 
decomposed into a product of three exponential opera- 
tors: 



ip{r,t + At) = exp 



d 2 
iAt^ 
2 



exp 



x exp 



d 2 
iAA 
2 

•a d 2 z 
\At^ 
2 



V(r,i). (12) 



Each exponential operator is approximated by Cayley's 
form as 

. , l + iAtd2/4 l + \Atd 2 /A 
^' t + A ^ l-iAt4/4 - l-iAt4/4 - 



1 + iAtd 2 /4 
l-iAt5 2 /4 



i>(r,t). (13) 



78iV floating operations are required to compute 
Eq. (|l^) ; where TV is the total number of grid points in 
the system. The norm and energy are conserved exactly. 
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By the way, a conventional method, Peaceman- 
Rachfold method Jl J^] , utilizes similar approximation ap- 
pearing on Eq. (J13 ) , which is a kind of the alternating 
direction implicit method (ADI method). However, by 
using exponential product, we have found that there is 
no need of ADI. This fact makes the programming code 
simpler and it runs faster. 



C. Static potential 

Next we consider a system subjected to a static exter- 
nal scalar field V(r). The TD-Schrodinger equation and 
its formal solution in this system are as follows: 



1 at 



-y + VSr) cir.t) 

A 



(14) 



ip(r, t + At) = exp [iAt— - iAtF(r)J ip(r, t) . (15) 

To cooperate with the potential in the framework of the 
formula described in the previous subsections, we have to 
separate the potential operator from the kinetic operator 
using Suzuki's exponential product theory P,|lO| as 



ip(r,t + At) 











exp 


[-4* 


exp 


iAt— 
2 J 



x exp 



i>(r,t) . (16) 



This decomposition is correct up to the second-order 
of At. The exponential of the potential is computed by 
just changing the phase of the wavefunction at each grid 
point. The exponential of the Laplacian is computed in 
the way described in the previous subsections. Each op- 
erator is exactly unitary, so the norm is conserved exactly. 
But due to the separation of the incommutable operators, 
the energy is not conserved exactly. Yet it oscillates near 
around its initial values and it never drifts monotonously. 
This algorithm is quite suitable for vector-type super- 
computers because all operations are independent by grid 
points, by rows, or by columns. The outline of this pro- 
cedure for a two-dimensional system is schematically de- 
scribed by Fig. [j]. 
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FIG. 1. The procedure for a two-dimensional closed sys- 
tem with a static potential. Here V shows the operation of 
the exponential of the potential, which changes the phase of 
the wavefunction at each grid point. K x and K y show the 
operation of Cayley's form along the x-axis and the y-axis re- 
spectively. They are computed independently by grid points, 
by rows, or by columns. 



The decomposition ( |l6| ) is a second-order one. Higher- 
order decompositions are derived using Suzuki's fractal 
decomposition || |Ojl3|. For instance, a fourth-order 
fractal decomposition S^At) is given by 

S 4 (At) = S 2 (sA<) S 2 (sA<) S 2 ((l - 4s)At) 

x S 2 (sAt) S 2 (sAt) (17) 

where 



(18) 
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D. Dynamic potential 

To discuss high-speed electron dynamics caused by a 
time-dependent external field V(r, t), we should take ac- 
count of the evolution of the potential itself in the TD- 
Schrodinger equation given as 

i^^=W(t)V(r,t) ; H(t) = -^ + V(r,t) . (19) 



The analytic solution of Eq. (|19|) is given by a Dyson's 
time ordering operator P as 



ip(r,t+ At) = Pcxp 



t+At, 



1 1 ^"{j-vM} 



1>(t, t) . 
(20) 



The theory of the decomposition of an exponential with 
time ordering was derived by Suzuki |fl2]| . The result is 
rather simple. For instance, the second-order decompo- 
sition is simply given by 



t + At) ~ exp 



x exp 



.At T „ At, 
-i—V(r,t+—) 





. Ai 


exp 


iAt— 

L 2 J 



-i^V(v,t+^) , •(!•./ 1 (21! 



and the fourth-order fractal decomposition is given by 

^(r, t + At) = S 2 (sAt; t + (l- s) At) 
x S 2 (sAt;t + (1 - 2s) At) 
x S 2 ((l-4s)At;t + 2sAt) 
x S 2 {sAt;t + s At) 
x S 2 (sAt;t) V(r,t) , (22) 



S 2 (At;t) =exp 



• At Tr/ At, 
-i—V(r,t+—) 

x exp 







exp 


iAt— 
2 J 



• At , At • 
-i—V(r,t+—) 



(23) 



These operators are also unitary. These procedures are 
quite similar to those of the static potential except that 
we take the dynamic potential at the specified time. 
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E. Periodic system 

In a crystal or periodic system, the wavefunctions must 
obey a periodic condition: 



^(r + R,i) = ^>(r, i) exp [i<? 



= k R 



(24) 



where k is the Bloch wave number and R is the unit 
vector of the lattice. The matrix form equation corre- 
sponding to Eq. (|^) in this system takes the following 
form: 



A -1 e +i< ^~ 
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Mt) 
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- Mt) - 



(25) 



These matrices have extra elements, so the equation can 
no longer be solve efficiently. 

We propose a trick to avoid this problem. We repre- 
sent the second spatial differential operator B 2 as a sum 
of two operators: 



This operation is exactly unitary and easy to compute. 

The exponential of 8 2 td is computed in the ordinary 
way. Thus the norm is conserved. We named d 2 d an 
"adhesive operator" because this operator plays the role 
of an adhesion to connect both edges of the system. The 
outline of the procedure for a two-dimensional periodic 
system is schematically described by Fig. g. 
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FIG. 2. The procedure for a two-dimensional periodic sys- 
tem. Here K x and K y show the operations of Cayley's form, 
and they operate as if this system is not periodic. X-adhesive 
and Y-adhesive mean the operations of the exponential of the 
adhesive operators along the x-axis and the y-axis, respec- 
tively. The operation of the adhesive operator needs only the 
values at the edges of the system. 



dt = d: 



xtd 



d 



x ml 



(26) 



F. Parallelization 



Multiplying by Ax 2 , the above representation reads in 
the matrix form: 



-2 1 

1 -2 1 

1 -2 

e +i< ^ 1 



-110 
1-210 
1-21 
1-1 



-10 





e +i0 0-1 



(27) 



The hrst matrix on the right-hand side, which corre- 
sponds to d 2 td , is tri-diagonal, and the second one, which 
corresponds to d 2 ad , is its remainder, and it has a quite 
simple form. The exponential of the second differential 
operator is decomposed by these terms: 



exp 



iAf 



dz 



exp 



iAt 



B 2 

u xa,d 



exp 



riAt 



B 2 



exp 



iAt 



xad 

(28) 



The exponential of d 2 ad is exactly calculated by the fol- 
lowing formula: 



exp 



-2iC 



-1 



(29) 



The adhesive operator plays another important role. 
It makes Cayley's form suitable for parallelization. Wc 
use the adhesive operator to represent the second finite 
difference matrix in the following way: 
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(30) 



The interior of the first matrix on the right-hand side 
is separated into two blocks, which means this system 
is separated into two physically independent areas. The 
second matrix, which is the adhesive operator, connects 
the two areas. A large system is separated into many 
small areas, and each area is managed by a single pro- 
cessor. Since the exponential of a block diagonal ma- 
trix is also a block diagonal matrix, each block is com- 
puted by a single processor independently. Data trans- 
mission is needed only to compute the adhesive oper- 
ator. The amount of data transmission is quite small, 
nearly negligible. The outline of the procedure for a two- 
dimensional closed system on two processors is schemat- 
ically described by Fig. |[ 
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FIG. 3. The procedure for a two-dimensional closed sys- 
tem on two processors. Adhesive shows the operation of the 
exponential of the adhesive operator for parallel computing. 
The operation of the adhesive operator needs only the values 
at the edges of the areas, so the data transmission between 
the processors is quite small. 



G. Adaptive mesh refinement 

It is necessary for real space computation to be 
equipped with an adaptive mesh refinement to reduce the 
computational cost or to improve the accuracy in some 
important regions. We improved the adhesive operator 
to manage a connection of between two regions whose 
mesh sizes are different, as illustrated in Fig. ^. 
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FIG. 4. An example of adaptive mesh refinement. The 
element in the left area is twice as large as that in the right 
area. The adhesive operator connects these areas. 

The second differential operator d 2 should be Hermite, 
but in this case the condition required for the matrix rep- 
resentation (dl)ij is given by 



(a 



for all 



(31) 



Considering this condition, an approximation of the 
second differential operator is given as 
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(32) 



The indices attached to this matrix indicate the corre- 
sponding mesh indices described in Fig. ^. This matrix 
is also divided into a block-diagonal one and an adhesive 
operator as 
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(33) 
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(34) 



The exponential of the adhesive operator is calculated 
using the following formula: 



exp 



iAt 

4Ax 2 



-1/4 1/8 1/8 
1/2 -1/2 
1/2 -1/2 

2ci —c\ —c\ 
— 4ci 2ci + c 2 2ci — c 2 
-4ci 2ci - c 2 2ci + c 2 



, (35) 



where 



ci 



C2 



exp 
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3i At 



V28A2; 2 
2i At 



V2 8Aa 
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1 
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(36) 
(37) 



In this way, it is found that the adhesive operator is 
important to simulate a larger or a more complicated 
system by the present method. 



III. APPLICATION 

In this section, we show some applications of our nu- 
merical method. Though these applications treat simple 
physical systems, they are sufficient for verifying the re- 
liability and efficiency of the method. Throughout this 
section, we use the atomic units (a.u.). 



A. Comparison with conventional methods 



As far as we know, the conventional methods of solv- 
ing the TD-Schrodinger equation are classified into three 
categories: 1) the multistep method [El, 2) the method 
developed by De Raedt [0 and 3) the method equipped 
with Cayley's form [||. 

In this section, we make brief comparisons between 
Cayley's form and other conventional methods by sim- 
ply simulating a Gaussian wave packet moving in a one- 
dimensional free system as illustrated in Fig. 0. 
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FIG. 5. The model system for comparison with con- 
ventional methods. 256 computational grid points are allo- 
cated in the physical length 8.0a.u. A Gaussian wave packet 
is placed in the system, whose initial average location x and 
momentum p are set at x — 2.0a.u. and p — 12.0a.u., re- 
spectively. 

The TD-Schrodinger equation of this system is simply 
given by 



(38) 



The wavefunction at the initial state is set as a Gaus- 
sian: 



tp(x,t = 0) = 



1 



v^nW 2 



exp 



\x - x 
AW 2 



IPoX 



(39) 



where W = 0.25a.u., x — 2.0a.u., p — 12.0a.u. 

The evolution of this Gaussian is analytically derived 

as 

1 

ip(x,t) 



^/2ttW 2 + (ir/2)(t/W) 2 

(X ~ X - p t) 2 

AW 2 + (t/W) 2 



x exp 



ip x 



(40) 



Therefore, the average location of the Gaussian (x(t)) is 
derived as if it is a classical particle: 



(x{t)) = (x(t = 0)) + Po t 



(41) 



This characteristic is useful to check the accuracy of the 
simulation. 

We use the second-order version of the multistep 
method and the De Raedt's method in order to compare 
with Cayley's form since Cayley's form is second-order 
accurate in space and time. 

The second-order multistep method we used in this 
system is given by 

d 2 

V>(i + At) = V(* - At) + i2A*-^ , (42) 
where d 2 is approximated by a finite difference matrix as 



(43) 



Extra memories are needed for the wavefunction at the 
previous time step ip(t — At). Though the time evolution 
of this method is not unitary, the norm of the wavefunc- 
tion is conserved with good accuracy on the condition 
that At/ Ax 2 < 0.5. This method needs only 107V float- 
ing operations per time step, which is the fastest method 
in conditionally stable methods. 
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Meanwhile, 
given by 



the second-order De Raedt's method is 



tp(t + At) = exp 



2 2 . 
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2 . 
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(44) 



where d 2 a and d 2 b are the parts of the second differen- 
tial operator and are approximated by finite difference 
matrices as below: 
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(45) 



(46) 



The exponentials of those matrices are exactly calcu- 
lated using the following formula: 

-2iC 



exp 



iC 



-1 1 
1 -1 



I 



1 



-1 1 
1 -1 



(47) 



The time evolution of this method is exactly unitary, and 
the norm is exactly conserved unconditionally. However, 
it seems that the accuracy tends to break down on the 
condition that At/ Ax 2 > 1.0. This method needs 187V 
floating operations per time step, which is the fastest 
method in unconditionally norm-conserving methods. 
Cayley's form with the finite difference method is given 

by 



i){t + At) = 



1 + iAt/4 d 2 x 



1 



(48) 



iAt/4 d 2 

where the spatial differential operator is approximated 
by the ordinary way in Eq. (f43|). 

The time evolution of this method is exactly unitary, 
and the norm is exactly conserved unconditionally. More- 
over, this method maintains good accuracy even under 
the condition that At/Ax 2 > 1.0. This method needs 
26./V floating operations per time step, which is the fastest 
method in unconditionally stable methods. 

We have simulated the motion of the Gaussian by those 
methods. First we show a comparison of Cayley's form 
with the conventional methods in the framework of the 
FDM. Figure |(] shows the time evolution of the error in 
the energy, which is evaluated by the finite difference 
method as described below 

e(t) = E(t) - E{t = 0) (49) 
1 Ar_1 

m = ~9AT Re £ #(<)(^-l(*) " 2 Mt) + ■ 

(50) 



i=0 



G 



The initial energy is evaluated as 73. 03a. u., though it is 
theoretically expected to be 74a. u. The ratio At/ Ax 2 is 
set at 0.5 to meet the stable condition required for the 
multistcp method. 

The energies violently oscillate in the results of the 
multistep method and De Raedt's method, as a result 
of the fact that these time evolution operators do not 
commute with the Hamiltonian. These energies seem to 
converge after the wave packet is delocalized in a uniform 
way over the system. Meanwhile, the energy is conserved 
exactly in the result of Cayley's form because Cayley's 
form commutes with the spatial second differential oper- 
ator which is the Hamiltonian itself in this system. 

Figure shows the relation of the time slice At to the 
error in the average momentum of the Gaussian, which 
is evaluated by the finite difference method as described 
below: 



e(M/A x *) = W = T))-( x (t = 0)) 



N-l 



- (p(t = o)) , 

(51) 



(x(t)) = AxJ^^M 2 (52) 

N-l 

(p(t)) = -Im]T , (53) 



Multistep 
De Raedt 
Cayley 



where T is a time span set at 0.4a.u. The initial momen- 
tum (p(t = 0)) is calculated as 11.7a.u., which is different 
from the theoretical value p — 12.0a.u. due to the finite 
difference method. 
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FIG. 6. Time variances in the energies computed by the 
three methods. The time slice is set at At = l/2048a.u. and 
the spatial slice is set at Ax = l/32a.u. so that the ratio 
At /Ax 2 is equal to 0.5. The energies violently oscillates in 
the result of the multistep method and De Raedt's method. 
Meanwhile, the energy is conserved exactly in the result of 
Cayley's form. 
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FIG. 7. Errors in the average momentum computed by the 
three methods in several time slices. The multistep method 
cannot be performed when At/Ax 2 > 0.5. The error of 
De Raedt's method is too large when At/ Ax 2 > 1. The 
error of Cayley's form is rather small. The spatial slice is set 
at Ax = l/32a.u. 



In the multistep method, the computation cannot 
be performed due to a floating exception, if the ratio 
At/ Ax 2 exceeds 0.5. In Dc Raedt's method, the error be- 
comes too large to plot in this graph if the ratio At/ Ax 2 
exceeds 1.0. Meanwhile, in Cayley's form, the error is 
not so large even if the ratio At/ Ax 2 exceeds 1.0. 

In this way, Cayley's form is found rather stable. 
Therefore, we can use a longer time slice than those of the 
other methods. And this Cayley's form becomes suitable 
for three-dimensional systems, potentials, periodic condi- 
tions, adaptive mesh refinement, and parallelizations by 
our improvements in this paper. 



B. Test of the adhesive operator 



To verify the reliability and efficiency of the adhesive 
operator for periodic condition and parallelization, we 
have simulated the motion of a Gaussian wave packet in 
a two-dimensional free system. As illustrated in Fig. ||, 
this system has periodic conditions along both the x-axis 
and the y-axis, and it is divided into nine areas, each of 
them is managed by a single processing element; the ad- 
hesive operator connects them. The initial wavefunction 
is set as a Gaussian given as 



1 



V2nW 2 



exp 



AW 2 



iPo ■ r 



(54) 



where r Q is set as the center of this system and p = 
(la.u., la.u.), W — la.u. The energy of this Gaussian is 
theoretically derived as 1. 0625a. u. 



7 





1 
















1 1 








I 


•E 


9 






I 


•E 


7 






I 


>E 


8 






I 


'E 


9 






I 


•E 


7 


— 








































































































I 


•E 


3 






I 


>E 


1 






I 


"E 


2 






F 


•E 


3 






I 


•E 


1 


— 








































































































I 


"E 


6 






I 


•E 


4 






I 


"E 


5 






I 


'E 


6 






I 


•E 


4 


— 








































































































I 


"E 


9 






I 


>E 


7 






I 


"E 


8 






I 


'E 


9 






I 


•E 


7 










































































































I 


•E 


3 






I 


'E 


1 






I 


•E 


2 






I 


•E 


3 






I 


•E 


1 






1 
















1 1 







FIG. 8. The model system for the test of the adhesive 
operator for periodic conditions and parallelization. This sys- 
tem is periodically connected and is divided into nine areas. 
Each area is managed by a single processing element. 32 x 32 
computational grid points are allocated in each area whose 
physical size is set at 8.0a.u. x 8.0a.u. The time slice is set at 
At = l/16a.u. 

Figure [)] shows snapshots of the time evolution of the 
Gaussian, which is observed to go through these areas 
smoothly. Figure ^ shows the evolution of the energy, 
which is observed to oscillate around its initial value. 




FIG. 9. Evolution of the density. The Gaussian is ob- 
served to go through these areas smoothly. 
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FIG. 10. Time variance in the energy. The initial energy 
is theoretically derived as 1.0625a.u., but it is evaluated as 
1.0553a.u. by the FDM. The energy oscillates near its initial 
value but never drifts monotonously. 

Second, we allocate 64 x 64 grid points only in the cen- 
tral area as illustrated in Fig. |ll|. We utilize the adhesive 
operator for the adaptive mesh refinement. Figure [l^ 
shows the snapshots, with the Gaussian going through 
these areas smoothly. Figure |l^ shows the evolution of 
the energy, which is observed to oscillate near its initial 
value. In this way, the reliability of the adhesive operator 
is proved. 
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FIG. 11. The model system for the test of the adhesive op- 
erator for the adaptive mesh refinement. This system is also 
periodically connected and is divided into nine areas. Each 
area is managed by a single processing element. The size 
of each area is set at 8.0a.u. x 8.0a.u. 32 x 32 computational 
grid points are allocated in each areas except the central area. 
The central area has 64 x 64 computational grid points, which 
makes it twice as fine as those of the other areas. The time 
slice is set at At = l/16a.u. 




FIG. 12. Evolution of the density. The Gaussian is ob- 
served to go through these areas smoothly. 
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1.063 




Time [a.u.] 

FIG. 13. Time variance in the energy. The initial energy 
is theoretically derived as 1.0625a.u., but it is evaluated as 
1.0591a.u. by the FDM. The energy oscillates near its initial 
value but it never drifts monotonously. 

C. Excitation of a hydrogen 

As the last application of the present method, we 
demonstrate its validity and efficiency in describing the 
process of photon-induced electron excitation in a hy- 
drogen atom in a strong laser field. The laser is treated 
as a classically oscillating electric force polarized in the 
z-direction: 



and the electric force. We follow the evolution for 32k 
iteration. 

Figure [l4| shows the time variance in the polarization of 
the electron. The oscillation of the polarization generates 
another electric field, which corresponds to a non-linearly 
scattered light from the atom. By Fourier-transforming 
the polarization along the time axis, we obtained the 
spectrum of the scattered light shown in Fig. [ll| 




Time [fs] 



FIG. 14. Time variance in the polarization of the electron. 
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E z = E a sin Lot 



(55) 



The spatial variation of the electric field of the light is 
neglected, because the electron system is much smaller 
than the order of the wave length. Then the interaction 
term of the Hamiltonian is approximated as 



Hint = ~eE z z . 



(56) 



In other words, we only take into account the electro- 
dipole interaction of the electron with the light, and ne- 
glect the electro-quadrapole, the magnetic-dipole, and 
other higher interactions. 

The amplitude E a is set at l/64a.u. = 0.80V/A, which 
is as strong as a usual pulse laser. The angular frequency 
to is set at 0. 3125a. u. = 8.5eV, less than the transition en- 
ergy between IS and 2P. Ordinarily, such low energetic 
electric force has no effect on the electronic excitation. 
But with such a strong amplitude, various nonlinear op- 
tical effects arc caused by the electron dynamics. 

We allocate 128 3 grid points in a 32 3 a.u. 3 cubic closed 
system. The hydrogen nucleus is located at the center of 
the system, and the nucleus potential is constructed by 
solving the Poisson equation in the discretized space to 
avoid the singularity of the nucleus potential. The 1S- 
orbital is assumed as the initial state of the wavefunction. 
Then we turn on the electric field and start the simula- 
tion. The time slice is set at 0.0785a.u. = 2.0 x 10~ 3 fs 
so as to follow the rapid variation of the wavefunction 



llli lli-lli. .,,!„. nll.illilljl mill I iIliI.Ii Jill, uJLl, ..ll 

5 10 15 20 25 



Photon energy [eV] 

FIG. 15. Spectrum of the scattered light generated by the 
oscillation of the electron. 

Several sharp peaks are found, which are interpreted as 
follows: The peak at 8.5eV comes from Rayleigh scatter- 
ing, whose frequency is identical with the injected light: 
lo. The peak at 10.2eV comes from Lyman a emission, 
which is generated by the electron transition from the 2P- 
orbital to the lS-orbital: LOL a - On the other hand, the 
peak at 12.1eV comes from Lyman (3 emission, which is 
generated by the electron transition from the 3P-orbital 
to the lS-orbital: lul i} ■ The peak at 6.8eV comes from hy- 
per Raman scattering, whose frequency is identical with 
2ia — loLc- Moreover the peak at 25.5eV comes from the 
third harmonic generation, whose frequency is identical 
with 3cj. 







The simulation is also performed for a different laser 
frequency; the injecting photon energy w is set at 10.2eV, 
which is the same as the transition energy between IS and 
2P. In this case the electron starting from a IS orbital is 
expected to excite to a 2Pz orbital. Figure |l^ shows the 
snapshots of the density during the simulation time span. 



One could obtain such behavior analytically by using 
perturbation theory; however, with the present method, 
we could directly calculate them without perturbation 
theory and without information on the excited states of 
the system. 



FIG. 16. Evolution of the density of the electron in the hy- 
drogen atom. The density starting from a IS orbital oscillates 
with time and becomes a 2Pz orbital. 

Figure [jj] and Fig. [l8] show the polarization and the 
spectrum, respectively. Three peaks are found, at 9.9eV, 
10.2eV, and 10.5eV. These peaks are derived from the 
theory of the Dressed atom or the AC stark effect as 
below: 



e£ (2P z |z|lS>, u, cj + e£ (2P z |z|LS) 



(57) 
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FIG. 17. Time variance in the polarization of the electron. 
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FIG. 18. Spectrum of the scattered light generated by the 
oscillation of the electron. 



IV. CONCLUSION 

We have formulated a new method for solving the 
time-dependent Schrodinger equation numerically in real 
space. We have found that by using Cayley's form and 
Suzuki's fractal decomposition, the simulation can be 
fast, stable, accurate, and suitable for vector-type su- 
percomputers. We have proposed the adhesive operator 
to make Cayley's form suitable for periodic systems and 
parallclization and adaptive mesh refinement. 

These techniques will also be useful for the time- 
dependent Kohn Sham equation, which is our future 
work. 
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